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Abstract 

This  paper  examines  an  inverse  problem  in  thermal  imaging,  that  of  recovering  a  void  in  a 
material  from  its  surface  temperature  response  to  external  heating.  Uniqueness  and  contin¬ 
uous  dependence  results  for  the  inverse  problem  are  demonstrated  and  a  numerical  method 
for  its  solution  developed.  This  method  is  based  on  an  optimization  approach,  coupled  with 
a  boundary  integral  equation  formulation  of  the  forward  heat  conduction  ])roblem.  Some 
convergence  results  for  the  method  are  proved  and  several  examples  are  i)resented  using 
computationally  generated  data. 
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1  Introduction 


Thermal  imaging  is  a  technique  of  recent  interest  for  the  nondestructive  evaluation  of  mate¬ 
rials.  This  method  attempts  to  characterize  the  i;iterna)  structure  of  a  sample  (perhaps  to 
locate  flaws-disbonds,  bubbles,  corrosion,  etc.)  by  using  its  surface  temperature  response  to 
an  external  heating.  Some  recent  work  and  techniques  in  this  subject  are  detailed  in  [3],  [4]. 
[5],  [7]  and  [9]. 

In  this  paper  the  problem  of  detecting  and  identifying  an  unknown  internal  void  in  a 
planar  domain  using  thermal  imaging  is  examined.  The  void  could  represent  a  defect  in  the 
material,  or  it  could  be  a  feature  which  is  supposed  to  be  present,  e.g.,  a  conduit,  whose 
location  or  geometry  is  to  be  assessed.  The  focus  is  on  the  case  in  which  the  thermal 
stimulus,  an  applied  heat  flux  at  the  boundary  of  the  sample,  is  periodic.  After  separating 
the  temporal  and  spatial  variables,  one  obtains  an  inverse  or  domain  identification  problem 
for  an  elliptic  equation.  Results  concerning  the  uniqueness  and  continuous  dependence  for 
the  inverse  problem  will  be  examined  and  an  algorithm  for  the  numerical  recovery  of  the 
void  will  be  presented.  This  algorithm  will  be  applied  to  examples  using  computationally 
generated  data  with  and  without  noise. 

The  outline  of  the  paper  is  as  follows.  Section  2  concerns  the  mathematical  formulation 
of  the  forward  heat  conduction  problem  with  periodic  heating  and  demonstrates  how  this 
leads  to  an  inverse  problem  for  an  elliptic  equation.  Section  3  contains  an  identification 
result  for  the  inverse  problem  and  shows  that  an  internal  void  is  determined  by  the  bound¬ 
ary  temperature  response  to  external  heating.  This  section  also  contains  results  concerning 
sensitivity  or  continuous  dependence  for  estimates  of  the  internal  void  based  on  the  bound¬ 
ary  measurements.  These  results  rely  on  reformulating  the  heat  conduction  problem  as  an 
integral  equation  on  the  boundary  of  the  sample.  Section  4  examines  a  numerical  method 
for  the  solution  of  the  forward  problem  based  on  the  boundary  integral  formulation  and  its 
incorporation  into  a  least-squares  routine  for  solution  of  the  inverse  problem.  It  is  shown 
that  with  reasonable  hypotheses  on  the  class  of  voids,  the  numerical  method  will  converge 
to  the  solution  of  the  least-squares  formulation  of  the  inverse  problem.  Section  5  presents 
the  results  of  this  algorithm  applied  to  computationally  generated  data. 


2  Mathematical  Formulation 


The  sample  (without  void)  to  l)e  tested  will  he  dt'iioted  1)V  ih  a  lauinded  region  in  IH^  with 
boundary.  The  internal  void  will  l)e  denoted  by  I).  wlitut'  1)  CC  witli  ('^  boundary. 
The  function  T{t..r)  will  denote  the  solution  to  the  heat  etjuation 
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where  ly  is  the  outward  unit  normal  vector  field  on  the  boundary  of  if  \  D.  k  is  tlu'  tlu'rmal 
diffusivity  of  if.  o  is  the  thermal  conductivity  of  if.  To  denotes  the  initial  temijerature  of 
the  region  and  <j  is  the  heat  flux  at  the  boundary.  Both  s  and  o  are  assumed  to  Ix'  known 
constants.  Of  course  it  is  assumed  that  g  is  not  identically  zero. 

The  heat  flux  g{.r.1)  will  be  assumed  to  be  periodic  in  tinu'  with  known  fnHjuency  ~  so 
that 

g{.rj)  =  Re{c‘-\j{.r)} 


for  some  complex  valued  function  </(.r).  Actually,  g  would  also  generally  include  a  steady- 
state  term  ar  well,  fuit  since  only  the  periodic  res]K)nse  is  of  intiuest.  this  term  can  lie  ignored. 
Under  this  assumption  one  can  sejiarate  variables  to  find  that  7'(.r./)  =  Re{7’(.r)f when’ 
satisfies 
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at  least  for  time  large  enough  so  that  the  initial  conditions  do  not  matti-r.  .Wjte  that  tlu' 
function  /’(.r)  solving  (2.1)  will  be  complex- valued,  consisting  of  a  n'al.  or  in  phasi’.  and 
imaginary,  or  out  of  phase  part. 


The  ii)verse  problem  of  iiitc'rest  is  the  foilovvitig:  (Jiv<‘fi  a  k;i(nvn  applieil  heal  Dux  //.  can 
the  shape  and  location  of  tlu'  void  D  be  nnicinelv  dt'tennined  from  nieasuri'inents  of  7  on 
the  boundary  of  Provided  a  uni(iueness  result  hohls.  on<'  would  also  lik('  to  know  wlu'tlu'r 
D  depends  contiinuuisly  on  nnuisurements  of  T  on  ()ih  that  is.  how  sensitive  estimates  of 
D  are  to  noise  in  the  data.  Finally.  <in<'  would  like  an  eflicient  computational  algorithm  lor 
recovering  an  estimate  of  D  from  actual  <lata. 

3  The  Inverse  Problem 

3.1  Uniqueness 

A  uniqueness  result  for  the  inverse  problem  l(dlows  (uisily  from  basic  facts  al)out  ('Hijrtic 
operators. 

Theorem  3.1  ( Ihriquoics!^)  Ld  D\  and  Di  ht  two  .'Oibdnmnins  of  U  and  T\  and  Ti  thf 
corresponding  solutions  to  equal  ion  (2.1)  with  nonzero  .Xeuinann  data  g.  Let  S  be  a  portion 
of  dXl  with  positire  raeasurc.  Then  T\  =7-2  on  S  implies  that  D\  =  Di. 

Proof:  The  functions  1\  and  Ti  have  the  .same  Neuma?m  data  g  and.  sincf'  7’i  and  T2 
agree  on  S.  the  same  Cauchy  data.  Fnique  continuation  tor  (  llijrtic  operators  imi)li('s  that 
T\  and  7-2  agree  on  UD2).  Then,  for  example,  the  function  T2  has  a  vanishing  normal 

derivative  on  the  region  D\  \  D2.  hence  is  constant  on  D]  \  1)2-  If  0\  \  D2  0  then  by  the 
maximum  principle  T2  must  be  constant  throughout  TI  \  Dz-  contradicting  g  /  0.  .\  similar 
argument  shows  that  if  772  \  Dj  0  then  T\  would  be  constant  on  Q  \  77].  again  (X)ntraflicting 
g  ^  0,  hence  /7|  =  D2. 

3.2  Boundary  Integral  Formulation 

In  order  to  investigate  continuous  dependemx'  and  numerical  nu't  boils  for  tiu'  rc'covi'ry  ol 
D  it  will  convenient  to  reforiTUilate  the  h<‘at  conduction  problem  as  a  lioundarv  integral 
equation.  This  offers  tlu'  advantage  that  one  oidy  has  to  solvi'  for  the  temi)erature  on  thi' 


boundary  of  the  sample,  rather  than  the  interior.  Since  the  boundary  is  the  only  place  the 
temperature  is  measured,  this  is  the  only  place  its  value  is  needed. 

Use  r(7-)  to  denote  a  fundamental  solution  or  CIreen’s  function  for  the  operator  (A  —  ^). 
Such  a  fundamental  solution  is  given  by 

r(r)  = 

=  -  — (kcr(ryuj//c)  +  ikc\{r  \J lo  j  k)) 

where  Kq  is  the  zero  order  modified  Bessel  function  of  the  second  kind  and  ker  and  kei  denote 
the  kelvin  functions.  Efficient  routines  for  computing  the  kelvin  functions  can  be  found  in 
[1].  Define 

r(x,?/)  =  r(|x  -  ?/j). 

The  function  T  satisfies  the  heat  equation  in  the  y  variable  for  fixed  x,  except  at  x  =  y, 
where  it  has  a  logarithmic  singularity.  Standard  potential  theory  arguments  (see  [6],  chapter 
3)  show  that  the  elliptic  problem  given  by  equation  (2.1)  can  be  formulated  as  a  boundary 
integral  equation, 

-^7’(x)+  [  T{y)d^^r{x,y)dSy=  I  Tix,y)g{y)  dSy  (3.1) 

2  J8(n\D)  Ja{Q\D) 

for  each  x  €  ^(D\  D)  where  is  the  normal  derivative  in  the  y  variable  and  dSy  is  surface 
measure.  We  will  use  K{x,y)  to  denote  the  kernel  d^^r{x,y)  for  x,y  E  d{il  \  D)  and  use  S 
to  denote  the  operator 

(S4>){x)=[  I<{x,y)<j)(y)dSy.  (3.2) 

Jd[Q\D) 

The  operator  S  is  bounded  and  compact  on  C{d{il  \  D)),  the  space  of  continuous  functions 
on  d{Q.\  D),  hence  —  i  /  +  .S'  is  a  second  kind  Fredholm  operator.  Uniqueness  of  the  solution 
to  equation  (3.1)  follows  from  uniqueness  of  the  solution  to  the  forward  problem.  By  the 
Fredholm  alternative,  the  equation  (  — ^7  +  S)(i>  =  g  is  solvable,  at  least  for  smooth  enough 
g.  In  particular,  (  — |  +  ‘S)“’  exists  and  is  bounded  on  C{d{il\  D)).  The  solution  to  equation 
(3.1)  yields  the  temperature  T{x)  on  80,  and  80.  If  needed,  the  temperature  for  x  G  U  \  7) 
can  be  found  from  the  relation  (Green’s  third  identity) 

T{x)  =  [  {T{y)8^^rix^y)-r{x^y)g{y))dSy. 

Ja{u\D) 
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3.3  Restrictions  on  the  Domain 

In  order  to  obtain  continuous  dependence  results,  a  few  restrictions  on  the  class  of  voids 
and  their  parameterization  are  needed.  First,  let  us  use  C^[0, 1]  to  denote  the  space  of 
functions  on  [0,  1]  where  the  endpoints  0  and  1  are  identified  with  each  other.  This  space 
can  be  normed  by  ||(?i||  =  supjgjo.i]  I/^I  —  2-  assume  that  D  depends  on  finitely 

many  parameters,  D  =  D{q)  with  q  ^  Q  CC  1R”‘  and  where: 

(a)  q\  =  (j(2  implies  D[q\)  —  D{q2)  (unique  parameterization). 

(b)  D{q)  C  ft'  CC  ft  ioT  q  ^  Q  (D(q)  stays  away  from  dfl). 

(c)  The  closed  curves  dD{q)  are  parameterized  as  x{q,t)  =  {x\{q,t),X2{q,t))  for  q  &  Q, 
0  <  t  <  1,  with  Xi{q,t)  a  function  of  t  for  each  q  £  Q  and  ^ 

bounded  away  from  zero.  Also,  the  map  q  — >  x{q,t)  is  continuous  from  IR”^  to  C^[0,  1]. 
Based  on  the  above  assumptions  it  is  not  difficult  to  show  that: 

(d)  For  each  x  €  dD{q)  there  exists  an  open  ball,  B^{x),  of  radius  e  around  x  with 
e  independent  of  x  and  q  such  that  the  curve  dD{q)  fl  B^{x)  is  parameterized  by 
(xi(<7,  t),  X2(<7,  t))  with  t  in  some  connected  interval. 

The  continuity  of  7  — >  and  compactness  of  Q  imply  that  the  norm  of  x{q^t)  as  a 

function  of  t  is  uniformly  bounded  over  q  ^  Q.  Also,  for  <7  E  Q  the  families  of  functions 
x{q,t),  x'{q,t)  and  x''{q,t)  are  equicontinuous  in  f,  where  the  prime  denotes  differentiation 
with  respect  to  t. 

Lemma  3.1  The  family  of  functions  {l\{x{q,s),x{q,t));  q  E  Q} ,  are  (unifomnly)  equicon¬ 
tinuous  in  s  and  t,  that  is,  for  each  e  >  0  there  is  a  S  >  0  so  that 

\I<{x{q,s),x{qJ))-  K{x{q,s'),x{q,t'))\  <  e 

for  all  (s,<)  and  (s',T)  with 

|.s  ~  s'j  <6,  \t  —  t'\  <  8, 


and  8  does  not  depend  on  s,  t  or  q. 


Proof;  For  brevity  let  us  suppress  the  depeudence  of  ./•()  and  /\’()  on  (j  and  write  simply 
for  K{x{s),  x{t)).  Also,  sinre  r(.r,/y)  =  log(j:r  -  //|)  +  (!{x.!;]  where  (!  is  smooth  in 
X  and  y,  we  will  prove  the  lemma  assuming  that  V{x,y)  =  log(|.r  —  //|);  it  is  easy  to  rheck 
that  the  smooth  term  makes  no  difference  in  the  proof. 

The  stated  regularity  for  l\  holds  on  the  compact  .set  {(s./)  G  [0,  1]  x  [(),  1],  |.s-/|  >., 
ty  G  f?},  where  f  is  any  number  greater  than  zero,  for  on  this  set  K  is  at  least  We  ik'cmI 
only  to  show  that  A’(,s./)  is  uniformly  continuous  on  th<'  st't  (.s, /)  G  [0.  1]  x  [0.  l].i.s-/|  <  (. 

Suppose  the  boundary  of  D  near  a  point  .ro  is  parameterized  by  ./■(/)  =  (.cif/).  .r^l/))  By 
making  an  appropriate  translation  and  rotation  it  tnay  lx*  assumed  that  ./'o  =  .;■(())  =  (0.0) 
and  that  the  boundary  is  oriented  so  that  the  unit  normal  outward  vector  in  (.ri..C2)  coor¬ 
dinates  is  (0,  1).  In  this  rase  Taylor’s  theorem  can  be  used  to  expand  x{t)  as 

.ri(/.)  =  (ifit  + +  R-\{i) 

•'■■2(0  =  + /^2(0 

where 


(iQ  =  ,ri(0) 

b  =  4(0) 

4(0  =  ^(4(O-4(0))/O  j  =  \.2. 

and  c  is  some  point  between  0  and  t.  The  functions  /?]  and  are  funct  ions  whose  norms 
can  be  bounded  in  terms  of  the  norms  of  .C]  and  .Ci-  The  unit  normal  vectors  satisfy 

^  .4(/)  = />/ +  4(0 
'^2(0^  =  -.r'(0  =  -<'o  -  -  4(0- 

at 

The  kernel  l\{sj)  is  then  given  by 

1  p/  /  X  in)—  -  1  (•'•Q-O  -  ■'•i(0)'^i(0  +  (-'•2(-s)  -  ■f2(t))i'2{J))(IR 

r.,  Lr,(s)-x,{t)y  +  {xA-s)~x-AI)V  <ll 

1  (.r,(.s)-.r,(/))4(O-('^-2(-O-.'-2(/)).0(0)  pj.jx 

'In  (•'■]  (.'i)  —  -r  1  (  /  ))^  +  ~  ■>'2{  0  )“ 

ti 


Substituting  the  al)ove  expressit)iis  f(.)r  ,r  i .  Xj. -c',  aiul  .r'^  givi's  for  (he  iiuiiH'iato!’.  after  ^oine 
simplification, 

nutn  =  — ^  “h  ~  n  ''u  +  ^  ^  ^  [/f  i  ( ■■' )  ^  ^ ) 

~  ( )  ~  I ) )  [Ud  +  e  1  /  +  ( f )]  —  ( ■‘^*  ~  I"  )H\{t  ]. 

Tlie  denominator  becomes  2-  (imes 


„o(/  +//,(..)-  /y,(/)  '+  -  t^)  +  IW) 


Taylor’s  theorem  implies  that 

=  Rj{t)  +  (-^  -  O/V'//)  + 

for  some  point  c  Iretwa'en  .s  and  I.  Substituting  this  into  the  ('xpif'ssion  lor  the  numerator 
gives 

mmi  =  — - — ^ — -  [— U(|6  +  ( 6/  +  H y{t)) li ((c)  —  (c/o  +  <; i /  +  H\{l )) (i\  li'iit )  ~-  )] 

where  c  and  c  are  betwe('n  s  and  /.  Doing  tlu'  same  for  th<'  (h'uominator  shows 
(U  n  =  'iTrl-s  —  ty  |((/u  +  /^,{  /  ))^  +  ( +  f  ’,.s  +  (  ^/j 


where  (\  and  ('2  functions  v)f  and  /  whicli  can  be  boiimlecl  in  (('tins  of  the  (’’  noiin  of 
x{t).  Cancelling  the  common  (s  —  /)^  factor,  the  kernel  <an  l)e  written  as 

/V(s  t]  ~  ^  — U06  +  [ht  +  )  -  [(/,;  —  iijl  —  +  (ijH'yjt)  —  ^ 

■l/T  [ui,  +  li\  (  /  )]'  +  [^^2^  n]'^  +  f  I  +  f  2^ 

Since  R\(0)  —  0.  the  keriud  is  bounded  through  ^  =  /.  Ir)r  the  denoininattu'  is  bounded 

away  from  zero  for  .s  and  /  in  a  sufhciently  small  neighborhood  of  zero.  In  fact,  by  raylor’s 

theorem  R\{f)  =  R'\{^^)  +  R'{{c)t  =  R"i<’)l  for  some  c  b('twe('n  zi'ro  ami  /.  In  addition,  since- 

the  functions  C]  and  ('2  can  lee  boumled  in  terms  of  the  f norm  of  .;(/)  (which  is  itself 

bounded  uniforndy  over  Q).  tin-  (h-nominator  may  be  boumled  uniformly  away  from  zeim 

for  .s  and  I  in  a  neighliorlueod  of  0,  with  the  m’ighi«irhoo<l  ami  bound  imh-pemh-nt  <>l  ./’n 

and  f/.  The  families  of  functions  •/’'(/)  and  as  well  as  R'J/}  ami  /)'”(/).  are  uniformly 

erpiicontimious  for  (/  €  Q.  ami  hence  st)  aia-leoth  thenmm-rator  and  ilenominat or  of  eipiat  ion 


(3.4).  Since  the  deiioininator  of  A’(.s./)  is  hounded  away  fiaun  zero,  it  follows  that  the  iainily 
of  functions  {/\  (x(</,  .s),  .r((/. /));  </  €  is  uniformly  e(|uicont  inuous  in  s  and  /.  ■ 

By  parameterizing  (.^ /)(</)  with  <  /  <  1  as  above  and  parameterizing  with  ./  (/). 

1  <  /  <  2,  one  can  identify  the  houndary  off7(fJ\/))  with  the  fixed  intt'ival  [0.2).  A  solution 
T(x)  to  equation  (3.1)  can  lie  identified  with  a  function  7'(/)  on  10.2)  1)\'  /  (/)  =  /  (./(/)). 
For  a  function  (f)  defined  on  [0,2)  h‘t  Oi  and  Aj  denote  the  restriction  of  o  to  [0,  1)  and 
(1,2),  respectively.  We  will  work  in  the  s|)ace  of  functions  o  for  which  Oi  is  c(uitinuous  and 
extends  continuously  to  [0,  1]  with  Oi(0i  =  Oifl)  and  for  which  O2  is  c(uit inuous  and  (-xtends 
continuously  to  [1,2]  with  Oj)!)  =  02(2).  Ue  will  denote'  this  space  hy  t^'[().2]  and  norm  it 
with  the  supremum  norm.  The  solutions  /’(/)  to  (3.1)  lie'  in  this  space.  Ont-  can  also  identity 
the  operator  S  (a  function  of  q)  with  an  integral  operator  on  t'’[0.2].  L<'t  us  write  l\(q.s.f) 
instead  of  i\{x{q.s).j-(q,t))  and  th'hiu' 

.S'((/)0(.s)  =  ^  {\{q,sj)0(t)  —  (lt.  (3..')) 

The  same  argument  given  in  the  proof  of  Lemma  3.1  sliows  (hat  /\’((/..s./)  is  uniformly 
equicontinuous  for  .s  and  t  in  [1,2),  i.e..  for  .r(.s)  and  j-{t)  on  the  houndary  of  H  and  q  6  Q. 
(actually,  l\  here  is  independent  of  q).  For  .s  6  [0.  1)  and  /  G  [1,2)  (e)r  vice-v<'rsa)  /\'(f/..s./) 
is  in  s  and  t,  since  in  this  case  x{s)  G  OD,  x{t)  G  OQ  and  hy  assumption  tiu'  houndarie's 
do  not  intersect,  so  again  K{q,s,t)  is  uniformly  equicontinuous  over  q  G  Q.  The  kerind 
K  will  have  a  jump  discontinuity  at  .s  =  1  or  /  =  1.  since  there  x[t)  jumps  from  OD  to 
(7n.  In  summary,  the  family  of  functions  {A'(<7,  .s, /);  q  G  Q},  is  uniformly  ('quicontinuous  on 
[p,  p+  1)  X  [<7,(7+  1)  p,q  =  1,2,  with  simple  jump  discontinuities  at  .s  =  1  or  /  =  1. 

Another  fact  worth  noting  is  that  the  map  q  — >  K{q. .•>.!)  is  confitiuous  as  a  map  from 
IR”*  to  6^[0,  2j  X  (7[0, 2],  that  is,  for  any  f  >  0  there  is  a  S  so  that  [(/  —  r/j  <  implies 

sup  |A'((7,s,/)  —  A"((7,.s, /|  <  f. 

s.i€[0,2] 

This  follows  directly  from  equation  (3.4)  and  the  fart  that  q  Xj(q.t).  q  — >  x''{q.f). 
q  —y  R'j{t)  and  q  — >  /?"(/)  are  all  continuous  as  maps  from  IR'"  to  the  space  of  continu¬ 
ous  functions. 
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3.4  Continuous  Dependence 

Based  on  tlie  above  assnm[)tions  it  is  possible  to  provt-  a  \'ersitai  cjf  cunt iimous  de|jeiideii<'e. 


Theorem  3.2  Ltt  (y,i  bt  a  sequitirt  i)i  Q  ami  I  (>1^)  Ibt  covvt  spoudnuj  .'^olitlnni  to  tqiiatioii 
{2.1}  with  I)  —  [){q,^).  Sup])os(  /’(</, i)  — ‘  I(q')  ill  \  D)]  for  souk  (/*  (E  ■  llnti 

<ln  q‘- 

Prt)of:  1  lie  first  task  is  to  show  that  tlie  map  </  — +  S(<j)  is  eont iuiious.  I.et  Ni  =  V(i/)  denote 
the  boundary  integral  opi-rator  for  [)  =  l)(q)  (eonsidererl  as  an  operator  on  r[().  2])  and 
■'^2  =  -^(V  "b  hq),  vvlien*  q  (z  Q  and  hij  is  some  small  pint urbatioti  in  </.  As  remarki'd  abo\c. 

q  — >  l\(q.s,t)  is  contitmoiis.  It  follows  that  q  — >  .'-'{q)  is  also  eontinuons  as  a  map  from  HP" 

to  the  spare  of  operators  on  ('’[0,  2j.  for  if  |A  (</  +  hq.  s.t)  —  h{q.  s.  /  )|  <  t  for  all  s.  /  G  [0.  2] 
then 

IK-SoK.s)  -  (>V>)(-‘*)||  <  /  \l\(q  +  hq.s.t)  -  l\{q.s.i)\\o{t)\ilt 

Jo 

for  some  constant  ('.  that  is.  ||>'((/  -irhq)-  >'(</) II  <  (’(  for  hq  suflieient !\'  stiiall. 

The  next  step  is  to  show  that  q  —*  (f  —  >'(</))“’  is  contitiuous.  IIk'  operator  (  — ^/  +  Pi ) 
is  inv'ertible  and  {  —  r^l  +  A'i)  van  be  inv<>rte<l  with  a  .Neumann  series  as  follows.  First  note 
that 

(--^l  +  S2r'  =  (-:^/  +  .S -(.s 

=  {--,f  +  >',)-'(/  +  +  .S)(Pi  -  Nb))"'. 

Let  H  =  ~{  —  jl  +  Pi)(.N|  —  .Pi)).  (liven  any  i  >  0  one  ran  ehoo.st'  a  |(*'f/|  suffirii-utly  small  so 

that  \\f{\\  <  f  and,  for  f  <  1. 

(-^/  +  .S'i)-'  =  (-1/  -t- 

=  {-~l  +  +  R  +  R'  +  ■■■) 

so  that 


aii<l 


( 


I  Inis  t  lu'  iiia|»  (/  ^  (  -  7  /  r  ■"'  (  </ 1 )  '  is  cdiiI  iinioiis  ami  so  I  he  iiia|)  ly  -•  /  I  </ )  —  (  -  ;  r  >'(  ly  i  '  (/ 
is  coiit  iiiuoMs  Iruin  Hr  ‘  to  (  ’'ll.  'jl. 

1()  Corn  I  )lt't  t‘  till'  pniul.  sii|)|)t>s<'  tlial  tin'  scijiicnct'  ly,  docs  not  coiivciiic  to  <j' .  Since  (J 
Is  (oiiipact,  some  sul)se(|',ience  of  <y,  eoiiveit^es  to  some  (j  e  (J .  It  can  l>e  assumed  that  this 
is  simpK  the  seipience  ly.,.  llo\ve\<'!.  eoiitmmty  o|  <y  -•  liij)  means  Unit 
iiiiphiii!'  /  liy')  =  /  ((j)  <iiid  eont  radiet  mu  the  uiii<pieiiess  I  heoreni  d.l. 


4  Numerical  Methods 

4.1  Nystrom’s  Method 


riie  t'oHowiiiu  ciimpiit  at  ioiial  apiuoach  to  the  iii\'eis<'  prohh'in  is  hasc-d  on  a  h'ast -scpiares 
(ormulatioii.  litidinu,  tlie  model  \oid  patainet<'rs  uhieh  Ik's)  lit  tlu'  iiK'asiiied  data  l»y  means 
ol  an  optimization  metlioil.  One  drawhaek  to  tliis  appi'oach  is  the  need  to  repoati'diy  soK’e 
the  lieaf  eoiidiict  ion  |)rohlem  (2.1).  If  is  thus  ad\' itita,u,<'<nis  to  have  an  ellieient  method  lot' 
soKinu,  this  eipiation.  1  lie  hoimdary  integral  equation  ap|)ioach  is  sneh  a  method,  and  in 
this  sect  ion  ’.ve  examine  a  te(  hnique  lor  its  sol  ut  ion. 

Idle  l>oundar\’  integral  lormnlatiiui  (d.l)  lor  the  solution  /  to  (-(piation  (2.1)  can  he 


writ  ten 


;'-7'(.s)  4-  /"  /\(,s, /)/•(/)<//  =  ,</(s) 

2  Id 


wlii're  /  (.s)  means  /’|,r(.s))  for  the  parameterization  .r(.s)  ot  ()([}  \  P).  Surtace  iiK'asiire 


■IS 


has  been  included  in  the  kernel  A  .  I.et  I ,  and  u,',  j  =  1 .  •  ■  •  .  n  deiiot e  the  nodes  and  weight s 
ot  a  (|uadrature  rule  comergent  on  the  sjotec  (''[0.2].  so  that 


^  u,', /(/..)  ^  /‘/(/) 

ml  ■" 


■II 


ID 


if  f{t)  is  smooth  enough.  .Actually,  we  will  consider  a  family  of  <iuadiaMire  nih's.  indexed  hy 
n,  with 


lim  =  f  f{l)(lt 

Jo 

tor  continuous  _/.  e.g..  the  n-node  composite  trapezt)idal  ride.  We  will  also  assume  that  tin' 
quadrature  rule  converges  uniformly  over  any  .set  of  ('quicont inuous  functions  in  f'[0.2]. 
so  that  if  /  G  .F. 


<  ( 


for  n  >  .V(f).  A  (()  independent  of  /.  The  n-nod<'  trapezoidal  rule  is  an  I'xample  of  such  a 
family,  or  more  appropriately,  tlu'  trapezoidal  rule  aiiplied  separati'ly  to  each  interval  [0.  Ij. 
[1.2]. 

Nystrom's  method  consists  of  replacing  the  integral  in  equation  (4.1 )  with  the  quadrat uri' 
rule  to  obtain 


—  :j7"„(.s)  +  ^  A (.s, /j)u;jT„(/^)  —  <i(s).  (4.2) 

Now  let  .s  =  A,^2W  •  •  to  obtain  the  i>  x  n  lini'ar  .systiun 


-  -Tn(t,)  +  E  =  .'/(/,)•  ( l.:i) 

j=\ 

The  idea  is  that  7'n(/i)  ~  T{t,).  .As  shown  in  [2],  each  solution  to  eijuation  (4.2)  h'ads  to 
a  solution  to  equation  (4.4)  am'  moreover,  each  solution  to  eipiation  (4.4)  corri'sponds  to  a 
unique  solution  to  etjuation  (4.2)  with  which  it  agri'es  at  the  nodes  /].•••./„.  Eiiuation  (4.4) 
is  the  system  which  is  solved  numerically  although  (4.2)  is  the  eijuation  we  will  usi'  for  the 
error  analysis. 

Write  equations  (4.1)  and  (4.2)  as 

+  =  2/ 

(--f  +  .s,)'/;,  =  1/ 

where  S,^  is  the  operator  in  eijuation  (  '.2).  Note  that  S„  is  coiiqiact  sine*'  it  is  a  liniti'  rank 
ojrerator  on  f’[0,2].  Ki’call  that  S  and  ilepeml  on  the  jiaranu'ter  (/. 

Theorem  4.1  I\,  — >  T  to  f‘[().2]  U.s  n  — >  oc.  Iht  coort  rqt  net  is  noifonii  nrtr  (j  G  (J ■ 
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In  proving  this  theorem  the  following  result  will  be  useful.  It  can  be  found  in  [2],  section  3.0. 


Theorem  4.2  Let  X  he  a  Banach  space,  let  S  and  (A  —  be  bounded  linear  operators 
on  X ,  with  A  7!^  0.  Let  T  be  a  bounded  linear  operator  on  X  with  the  property  that  either  A 
is  an  eigenvalue  ofT  or{X  —  T)~^  exists  (ifT  is  compact,  this  is  satisfied).  Further,  assume 


\\{T-S)T\\< 

Then  (A  —  T)~^  exists  on  X  and 


A 


ii(A-.^)-‘ir 

l  +  ||(A-5)-' 


|a|-||(a -.?)->  Ilii(r-,s’)rir 

Furthermore,  if  (A  —  .S')/  =  5  and  (A  —  T)h  =  y  then 

I' -  ''lA|-||(A-.9)-’||||(r-.9)T||‘ 


Proof  of  Theorem  4.1  The  proof  of  Theorem  4.1  is  simply  an  application  of  Theorem  4.2 
with  X  =  C[0,2],  A  =  —1/2,  S  =  —.S'  and  T  =  — .9„.  In  the  previous  section  it  was  shown 
that  the  map  q  — >  (  —  5/  +  S{q))~^  is  continuous,  hence  ||  —  A/  +  'S(7))”'||  can  be  uniformly 
bounded  over  Q.  In  order  to  complete  the  proof  of  Theorem  4.1  it  must  be  shown  that 


||(.S'-.S'„).9,1|  0  (4.4) 

l|(■S'-.S'n)^/|(  ^  0  (4.5) 


as  n  — >  oc,  uniformly  for  q  ^  Q.  This,  in  conjunction  with  Theorem  4.2,  will  show  that 
II T  —  T„|j  — >  0  uniform!;  in  q. 

From  the  argument  given  in  section  3,  for  .s  6  [0,2]  and  q  ^  Q,  the  kernel  K{q,s,t)  is 
uniformly  equicontinuous  in  the  t  variable.  The  uniform  convergence  of  ||(.S'  —  .S'n)5||  to  zero 
follows  from 


(•S' -  .S'„)(7(.s)  =  f  {K{sJ.)g{t)dt 

j=i 

<  ^{n) 


with  f(n)  ^  0  as  n  — >•  00,  independently  of  q  and  .s.  Here  we  have  used  the  fact  that  A'  is 
equicontinuous  in  t  and  the  assumption  that  the  integration  rule  converges  uniformly  over 
equicontinuous  sets  in  C[0,2]. 
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The  case  for  the  convergence  of  ||(.S'  —  S'„).S'n||  is  similar.  Let  cf)  be  a  function  in  C''[0,2] 
with  ||</)||  =  1.  Since  K  is  uniformly  equicontinuous  in  .s  and  t, 

Sn(l>{s  +  As)  -  +  As,tj)  -  K{s, 

J=1 

where  f(As)  does  not  depend  on  s,  t  or  q.  For  convergent  integration  rules  the  sum  I2j-i 
is  bounded  in  n  (see  [2],  part  I,  section  4,  theorem  7),  so  that  Sn(f>  is  an  equicontinuous 
function,  independent  of  q.  The  rest  of  the  argument  is  as  in  the  previous  paragraph  but 
with  g  replaced  by  Sn(l>-  This  shows  that  ||(.S'  —  0  uniformly  for  q  £  Q  and 

completes  the  proof  of  Theorem  4.1. 

4.2  Application  to  Inverse  Problem 

Let  us  suppose  that  the  temperature  data  for  the  inverse  problem  consists  of  point  measure¬ 
ments  Ti  at  points  x,  on  the  boundary  of  fi,  i  =  1,  •  •  • ,  M.  A  reasonable  way  to  attempt  a 

recovery  of  the  unknown  region  D  is  to  define  the  functional 

M 

JM  =  ’EiTMin)  -  T.f 

i=l 

and  seek  an  estimate  of  D{q)  as  the  solution  to 
(IDP)  minimize  J{q)  for  q  £  Q 


In  practice  the  solution  to  this  problem  will  involve  not  the  true  temperature  T{q)^  but 
the  n-node  Nystrom  approximation  Tn{q).  The  actual  minimization  problem  solved  will  be 

(IDP)”  minimize  J"(9)  for  q  £  Q 


where 

M 

•n?)  =  Dr.('7)(ii)-r,)''.  (4.6) 

i=I 

Based  on  the  analysis  of  the  convergence  of  T„  to  T,  the  following  theorem  can  be  stated. 
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Theorem  4.3  Ld  q,^  hr  a  solution  to  (IDPy.  Thfn  as  n  — +  oc,  sonii  subsi  qiu  ixf  of  (/„ 
court  rgt  s  to  q*  G  Q-  Moirortr.  q"  is  a  solution  to  (ll)l’). 


Note  that  if  the  l)omKlary  data  7’,  are  ronsist(>nt  with  some  Diq')  for  tf  G  Q  then,  since' 
there  is  a  unique  such  c/*,  one  can  talk  about  the  unique  solution  to  (IPP)  and  a  suhsi'eiuence 
of  </„  will  converge  to  this  q‘ . 

Proof  of  Theorem  1.3:  Some  sidrsequenct'  of  q^  converges  te)  a  point  q'  G  Q  by  virtue'  e)f 
the  fact  that  Q  is  compact.  Let  r/,,,  be'  anv  .seH|uene  ('  in  Q  e  e)nverging  te)  e/”'.  From  the'  uniform 
convergence  of  7k  to  /'  over  Q  we  e  an  e'eiiiclude  that 


lim  ./"(e/m)  = 


For  any  q  ^  Q  we  have 


■r(q„)  <  •/"('/)• 

Taking  the  limit  over  u  and  using  (4.7)  shows  that 

JUn  <  Aq) 


for  all  q  G  Q,  i.e.,  e/'  solves  (IDP). 


(1.7) 


5  Implementation  and  Examples 

5.1  Introduction 

In  this  se’ction  the  recoveuy  algorithm  is  implemeuite’el  and  the  re'sidts  crl  some  nume'iical 
experiments  are  de’seribed.  One  moelilication  is  maele  to  the  assumptions  e)t  the  luevious 
sections;  for  experimental  work  it  is  more  conve'iiient  to  use  a  sample'  P  which  is  re'ctaugular. 
thus  the  boundary  of  (7  will  not  lie  This  make'  little  elilh're'nce'  tc;  the  inte'gral  ('epiation 
formulation,  for  the  oiierator  (/  —  S{q)),  while  no  k)nger  a  second  kind  Fre'dluelm  oi)('rator. 
is  in  fact  a  small  pe’rturleation  ot  such  an  ojx'rateer.  One'  can  still  e'stablish  the'  e'xist e'lie'e'  and 
boundeehiess  of  the  inverse  (/  —  S{q))~^:  se'e  [10]  for  me»re'  eh'tails.  The  re'st  eef  the'  argume'nts 
are'  unchanged.  The  solutieen  T  on  the-  benmelary  e)f  is  still  e  eeiit  inuous.  althoeigh  it  will  leave' 
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‘Vorners”  at  the  corresponding  corners  of  fi.  For  numerical  purposes  it  is  thus  Ijeneficial  to 
use  a  quadrature  rule  which  allocates  more  nodes  near  these  corners,  rather  than  >miformly 
over  the  boundary  of  ft.  Note  that  it  is  still  assumed  that  ()D  is  C^. 


One  other  modification  will  also  be  made:  The  heat  flux  ()  will  be  a  point  heating  source, 
so  that  ^(j-)  =  8p,  a  delta  function,  where  P  is  the  point  at  which  the  heat  is  applied.  While 
this  flux  g  does  not  live  in  the  space  of  functions  in  which  we  have  been  working,  as  will 
be  shown  one  can  analytically  remove  the  .singularity  from  the  problem  and  work  with  a 
smooth  remainder  term  which  fits  the  hypothe.ses  made  so  far.  In  the  examples  that  follow 
we  use  only  the  imaginary  or  out  of  phase  portion  of  the  periodic  temperature  response  T . 
This  part  of  the  temperature  is  continuous,  even  through  the  singularity  at  the  point  heating 
source.  The  real  or  in  phase  part  of  T,  however,  has  a  logarithmic  singularity  at  the  point 
source.  In  reality  one  does  not  have  a  point  source,  but  for  a  heat  source  concentrated  in 
a  small  region  the  in  phase  response  in  the  region  of  the  source  varies  radically  with  the 
heat  source  “footprint”,  which  causes  problems  in  a  reconstruction  algorithm.  In  practice 
therefore,  one  excludes  the  in  phase  data  near  the  source.  For  convenience,  we  simply  work 
only  with  the  out  of  phase  portion. 

One  final  remark;  The  linear  system  obtained  via  the  n-node  NystronTs  method  can  be 
written 

^(</)7’n(</)  = /(</)  (5.1) 


where  A((i)  is  the  matrix  generated  by  the  discretization  and  b{q)  is  the  right  hand  side  of 
the  discretized  integral  equation;  both  depend  on  the  parameter  g  which  defines  the  void. 
As  a  result,  the  temperature  Ty^  also  depends  on  g.  The  Levenberg-Marquardt  optimization 
routine  used  requires  the  derivatives  of  the  functional  J{g)  with  respect  to  g  which  in  turn 
requires  On  can  obtain  these  derivatives  by  differentiating  equation  (5.1)  with  respect 


to  g  to  obtain 


dj  dA 

itr = 5^  ■ 


(,5.2) 


Both  A{g)  and  f{g)  are  differentiable  with  respect  to  g  and  if  A{g)  is  not  singular  then  the 
above  computation  is  valid,  that  is.  exists  and  therefore  satisfievs  {'qiiation  (-5.2).  Once 


Tn  has  been  obtained,  can  be  obtained  from  (-5.2);  in  fact,  the  work  done  in  solving  the 
equation  (LIT  decomposition)  can  be  re-used  in  computing  the  derivati v('.  It  should  be  noted 


that  the  temperature  satisfies  equation  (5.2)  at  both  the  finite-dimensional  (matrix)  level 
and  infinite  dimensional  (integral  equation)  level,  i.e.,  ^  also  exists. 


5.2  Sample  Geometry  and  Implementation 

The  geometry  used  for  the  examples  is  illustrated  in  Figure  1.  The  lower  left  corner  of  the 
rectangular  sample  has  (xjyXj)  coordinates  (0,0).  We  will  examine  the  case  where  the  voids 
are  disks  of  unknown  center  and  radius,  D  =  D{q)  with  q  —  (a,  b,  r)  where  (a,  b)  is  the  center 
of  the  disk  in  (xi,X2)  coordinates  and  r  is  the  radius.  The  boundary  of  D  is  parameter¬ 
ized  by  ;ci(i)  =  a-|-7'cos(27ri),  X2(t)  =  6+r  sin(27rf),  0  <  t  <  I.  If  the  length  (rq  axis)  of  Q  is  £ 


Point  heat  source 


Figure  1:  Sample  geometry. 

and  the  height  {x2  axis)  is  H  then  the  boundary  of  $7  is  parameterized  from  t  =  \  to  t  =  2 

by 


/(^)  = 


(4££0) 

l<t<l 

4  -  ^  2 

(£-(/-§),//) 

2  <  <  <  I 

2  —  ‘  ^  4 

(0, // -  4//(<  -  D) 

With  appropriate  bounds  on  the  range  of  a,  b  and  r,  it  is  simple  to  check  that  this  cla.ss  of 
voids  and  parameterization  satisfies  the  properties  of  section  3.  The  dimensions  of  the  sample 
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for  the  present  example  will  be  1.27  cm  in  length  by  0.32  cm  in  height.  All  references  to 
coordinates  and  sample  geometry  will  be  in  centimeters.  The  thermal  parameters  correspond 
to  aluminum.  As  mentioned  above,  the  heat  flux  will  be  a  point  source  applied  at  a  points 
with  (xi,X2)  coordinates  (xj,  0.32  cm)  on  the  top  of  the  sample.  This  source  will  have 
unit  (one  watt)  power.  Note  that  for  the  full  time-dependent  problem  this  means  that  the 
variation  in  heat  flux  is  one  watt. 

The  quadrature  rule  used  for  Nystom’s  method  for  the  solution  of  equation  (3.1)  is  just 
a  version  of  the  midpoint  or  trapezoidal  rule  with  variably  spaced  nodes.  Specifically,  let 
j{x)  =  for  0  <  X  <  |.  The  variable  p  is  a  parameter  to  adjust  the  spacing  of  the 

nodes.  Allocate  N  nodes  f,  [N  even)  and  weights  Wt  on  the  interval  [0, 1]  by 

1  1  —  (N-1+1.  1  =  ^  +  1,...,^ 

The  corresponding  weights  are 

2(^1  +  ^2),  i  =  1 

~  ?  =  2, . . . ,  jV 

1  ~  +  f/v),  i  =  N. 

This  rule  allocates  nodes  and  weights  to  (0, 1),  excluding  the  endpoints.  For  p  =  1  it  is  just 
the  midpoint  rule  on  [0, 1].  Choosing  p  >  1  causes  the  nodes  to  “bunch  up”  near  0  and  1. 
For  each  of  the  four  sides  of  dfl  the  nodes  are  allocated  by  mapping  the  above  nodes  on 
[0, 1]  to  the  corresponding  interval  in  f,  e.g.,  the  nodes  are  allocated  on  the  bottom  of  ff  by 
transforming  the  nodes  on  [0,  1]  as  <  — >  |  1  with  the  corresponding  change  in  the  weights. 

The  parameter  p  is  chosen  to  be  2.0;  this  allocates  more  nodes  close  to  the  corners  of  fl. 
On  the  void  boundary  dD  the  nodes  are  evenly  spaced  with  weights  corresponding  to  the 
trapezoidal  rule  on  [0,  1].  The  typical  number  of  nodes  used  is  10  to  30  on  each  side  of  fi 
and  20  to  40  on  the  void  boundary. 

In  order  to  deal  with  the  singular  boundary  heating,  first  note  that  if  r(x,p)  is  a  funda¬ 
mental  solution  to  L  =  (A  —  ^)  then  L^r(x,  p)  =  0  for  p  ^  x,  where  Ly  means  the  operator 
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applied  in  the  y  variable  Moreover,  if  P  G  OQ  and  T{y)  =  —'>T[P.y)li\  then  it  can  be  shown 
that 

ad^^f{y)  =  6p  -  20^^  r  ( P,  y )  on  iKl . 

that  is,  for  delta  function  heating  T  is  the  solution  uj)  to  a  nonsingular  remainder  term. 
The  remainder  term  on  the  right  hand  side  is  actually  continuouN  on  PU,  or,  more  precis(>ly. 
extends  continuously  through  y  =  P.  Thus  one  can  analytically  remove  the  singidarity 
associated  with  the  point  heating  and  simply  .solve  for  the  .smoother  nmiainder  term,  ddie 
forward  heat  conduction  problem  for  the  example  is  then 

(A  -  —)T  =  0  in  n\D 

K 

ad^T  =  2d^^r{P,y)  on  (){Q\D). 

The  full  solution  with  delta  function  heat  flux  would  be  given  by  the  sum  l'(y)  +  T{y).  This 
solution  has  a  logarithmic  singularity  in  its  real  part  and  has  smooth  imaginary  part. 

For  reconstruction  purposes,  the  temperature  in  the  present  examples  will  only  be  mea¬ 
sured  along  the  top  of  the  sample,  so  that  the  least  squares  fit  functional  (4.6)  includes  only 
this  data.  The  heating  frequency  is  in  the  range  of  1  to  5  Hz.  For  each  heating  source  the 
temperature  response  is  measured  at  40  equispaced  points  along  the  top  of  the  sample. 

5.3  Strategy 

One  of  the  necessities  of  an  optimization  approach  is  that  one  have  a  reasonable  initial  guess 
at  the  true  void  before  beginning  the  optimization  procedure,  or  el.se  risk  become  trapped  in 
a  local  minimum  far  from  the  “true”  solution.  This  is  particularly  applicable  in  the  present 
case,  as  illustrated  by  the  following  figures.  The  sample  with  void  is  shown  in  Figure  2.  The 
sample  is  again  aluminum  with  the  same  dimensions  as  in  Figure  1.  The  true  void  (solid 
outline)  has  a  radius  of  0.06  cm  and  is  centered  at  (;r],  j’2)  coordinates  (0.88  cm,  0.24  cm). 
The  heat  source  is  applied  directly  on  top  of  D*  at  d  Hertz.  The  prospective  void  D  is  fixed 
to  have  the  same  radius  and  x-^  coordinate  as  D*\  its  .rj  coordinate  is  allowed  to  vary  from 
0.15  cm  to  1.15  cm.  The  value  of  the  functional  J{q)  as  a  function  of  .ri  is  shown  in  Figure 
3.  The  functional  is  of  course  zero  when  the  x\  coordinates  of  D  and  D*  coincide  and  the 


Heating 


1.27  cm 


Figure  2;  Sample  geometry  for  least-squares  fuurtioual  exami)le. 


Aesiauoi  vs.  void  x  — coo'dinote 


Figure  3:  Functional  J{q)  versus  D  X]  coordinate. 

functional  rises  steeply  as  one  moves  away  from  the  minimum.  However,  if  one  used  an 
optimization  method  with  an  initial  guess  which  was  far  from  the  correct  value  (.Ci  <  0.5 
cm)  then  the  optimization  routine  would  probably  not  be  successful,  for  in  this  region  the 
functional  is  almost  flat-actually,  it  slopes  slightly  away  from  the  minimum.  This  illustrates 
the  nf'ed  for  a  reasonable  initial  guess  at  the  x\  coordinate  of  the  void  D*. 

In  the  previous  example  the  heat  source  wa*^  applied  directly  on  top  of  the  true  void 
D* .  In  reality  if  a  single  heat  source  is  applied  it  will  not  likely  fall  on  top  of  or  even  near 
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D*.  Figure  4  illustrates  the  same  situation  as  Figures  2  and  3  hut  with  the  heat  sourre 
far  (xi  coordinate  0.35  cm)  from  D*.  Here  the  situation  is  even  worse,  for  now  the  least- 
squares  functional  has  many  local  minima  to  trap  any  optimization  method  started  with 
an  xi  coordinate  far  from  D*.  Also,  since  the  heat  source  does  not  “illuminate"  D*,  the 
functional  is  very  flat  near  the  minimum.  In  the  presence  of  noise  one  would  not  be  able  to 
locate  this  minimum  with  any  accuracy. 


Residual  vs.  void  x  — coordinate 


Figure  4:  Functional  J{q)  versus  D  xj  coordinate. 


This  leads  to  the  following  strategy  for  locating  a  void,  illustrated  in  Figure  5:  Apply 
the  heat  source  at  a  number  of  different  points  along  the  top  the  sample.  For  each  different 
heating  location,  take  the  corresponding  temperature  measurements  along  the  top  of  the 
sample.  As  one  passes  the  heating  source  over  a  void  it  will  be  detected  by  a  change  in  the 
temperature  response.  Suppose  this  occurs  when  the  heating  source  has  an  Xj  coordinate 
equal  to  a.  Then  begin  the  optimization  with  initial  guess  Xi  =  a  and  X2  and  r  anything 
reasonable.  This  should  provide  a  reasonable  approximation  to  the  correct  Xi  coordinate  of 
the  void;  the  initial  guess  at  X2  and  r  is  much  less  crucial. 
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Figure  5:  Strategy  for  applying  heat  sources. 


5.4  Results 

This  procedure  is  illustrated  by  the  graphs  in  Figure  6.  The  actual  void  D*  is  located  at 
(x),X2)  coordinates  (0.77  cm,  0.24  cm)  with  radius  0.06  cm.  A  total  of  9  point  heating 
sotirces  at  3  Hertz  were  applied  on  the  sample  top  surface,  equispaced,  spanning  one  half 
of  the  sample  length,  with  source  5  centered  on  the  top  surface.  The  graphs  show  the  out 
of  phase  or  imaginary  portion  of  the  temperature  response  for  each  source  location.  The 
temperature  response  changes  most  rapidly  and  peaks  between  heating  locations  3  and  4. 
Contrast  this  to  the  temperature  response  when  no  void  is  present,  shown  in  Figure  7;  the 
response  does  not  change,  since  nothing  under  the  heating  source  changes  as  the  source 
moves.  A  Levenberg-Marquardt  algorithm  implemented  along  the  lines  in  [8]  was  used  to 
recover  an  estimate  of  the  void  Z)*,  using  only  the  peak  temperature  response  data,  heating 
sources  3  and  4.  The  initial  guess  at  the  Xj  coordinate  was  halfway  between  these  sources 
at  0.75  cm.  The  12  and  radius  initial  guesses  were  0.16  cm  and  0.1  cm,  respectively.  The 
optimization  code  converges  to  the  correct  void  in  1 1  iterations,  reducing  the  residual  from 
0.154  to  less  than  0.002.  However,  no  noise  was  present. 

As  a  more  realistic  example,  we  take  the  same  geometry  and  void  as  the  previous  ex¬ 
ample  but  with  20  percent  zero-mean  gaussian  noise  (measure  as  percent  of  signal  RMS 
value)  added  to  the  temperature  response.  The  responses  are  shown  in  Figure  8.  The  largest 
response  is  for  position  3.  The  data  for  heating  positions  2,  3  and  4  were  then  used  in  the 
optimization  with  initial  guess  Xj  =  0.79  cm  (the  Xi  coordinate  for  source  3),  x-2  —  0.16  cm 


21 


and  r  =  0.1  cm.  The  initial  residual  is  0.406.  The  optimization  routint'  re(lnces  this  to  0..4Ni 
in  12  iterations.  The  final  estimate  for  the  void  is  .r,  =  O.To.j  cm.  ,('2  =  0.22  cm.  r  =  0.067 
cm.  The  true  void  is  shown  as  the  solid  outline  in  Figure  9  and  tlu'  ami  recovered  ('stimatc 
as  the  dotted  outline. 


Poii.  7  t-’05.  8 


Figure  6;  Out  of  phase  temperature  for  varying  heat  source  locations. 


Degrees  C  Degrees  C  Degrees 


F’os. 


O.IOi - 

0.00  - - 

-0.10-  '\x 

-0.20  - 

-0.30  — I — ' — i — I — ^ 

0.0  0.6  1.2 
Top  of  sample,  cm 

Pos.  4 

O.IOi - ^ - 

0.00  — ^  ^ 
-0.10  -  \y 
-0.20  - 

-0.30  I— J ^ — ' — ' — ^ 
0.0  0.6  1.2 
Top  of  somple,  cm 

Pos.  7 

0.10 
0.00 
-0.10 
-0.20 


-0.30 

0.0 


o  O.C 


5  -0.10  - 


rji 

a> 

rj  -0.20 


-0.30  >- 
0.0 


D.O  0.6  1.2 

Top  of  sample,  cm 

Pos.  5 


.OJO 
o  0.00 

in 

S  -0.10 

rji 

Q  -0.20 


Top  of  sample,  cm 


O  0,00 
in 

S  -0.10 


-0.30  ' — ' — ' — “ — ' — ' — ^ 
0.0  0.6  1.2 
Top  of  sample,  cm 


0.00  -- 

-0,10  - 

-0.20  - 

-0.30  — 

0.0 


Top  of  SCogjUi, 
i’C';.  f> 


0.10 

0.00 

-0.10 

-.0,20 


.0,30 1 — ^ ^ — I — 1 — I — D 
0.0  0.6  1.2 
Top  of  sample,  cm 

Pos.  8 


-0.30 

0.0 


Top  of  scmpie,  cm 
l-’O'S,  9 


0.10  r 
0.00  - 
-0.10  - 
-0.20  - 


-0.30'— J — 
0.0 

Top  oi 


Figure  7:  Out  of  phase  temperature  for  varying  heat  source  loea!  ions,  no  void 
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Figure  8:  Out  of  phase  temperature  for  varying  lieat  source  locations.  20  i^ercent  noi 


6  Conclusion 


The  problem  of  recovery  a  void  in  a  material  sample  based  on  the  sample’s  surface  tem¬ 
perature  response  to  external  heating  has  been  considered.  Uniqueness  and  continuous 
dependence  results  have  been  shown,  and  an  optimization  based  algorithm  for  recovering 
an  estimate  of  the  void  has  been  demonstrated.  Much  of  this  work  rests  on  reformulating 
the  heat  conduction  problem  as  a  boundary  integral  equation,  which  provides  a  means  of 
rapidly  solving  the  heat  conduction  equations.  The  algorithm  was  run  for  the  simple  case 
of  a  circular  void  and  computationally  generated  data.  This  algorithm  exhibits  some  of  the 
problems  that  optimization  appioaches  are  heir  to,  specifically,  the  likelihood  that  a  poor 
initial  guess  will  not  converge  to  a  global  minimum,  and  some  strategies  for  overcoming  this 
have  been  described. 

Collection  of  actual  data  from  a  experimental  setup  at  NASA  Langley  Research  Center 
has  already  been  performed.  Preliminary  analysis  of  this  data  shows  good  agreement  with 
the  heat  conduction  model  and  the  ability  to  actually  recover  subsurface  voids.  The  analysis 
of  this  experimental  data  will  be  reported  elsewhere.  Of  interest  for  future  research  is  a  study 
of  the  sensitivity  of  this  thermal  technique,  e.g.,  the  size  of  voids  which  can  be  detected  at 
various  depths.  Techniques  for  voids  of  different  shapes  (especially  cracks  or  disbonds)  should 
also  be  examined.  The  heat  conduction  model  could  also  be  improved;  zero  boundary  flux 
away  from  the  source  becomes  unrealistic  at  low  frequencies.  We  would  also  like  to  pursue  a 
full  three-dimensional  heat  conduction  model  leading  to  a  two-dimensional  boundary  integral 
formulation.  Such  a  formulation  would  require  a  finite  or  boundary  element  technique  for 
the  ’  neiical  solution,  rather  than  Nystrom’s  method. 
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